home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / lmitool / bigm.sci next >
Text File  |  1999-09-16  |  3KB  |  105 lines
  1. function [x,Z,z,ul,iters]=bigM(F,blck_szs,c,x0,M,nu,abstol,reltol,tv,maxiters);
  2. // [x,Z,z,ul,iters]=bigM(F,blck_szs,c,x0,M,nu,abstol,reltol,tv,maxiters);
  3. //
  4. // minimize    c^T x 
  5. // subject to  F(x) = F0 + x1*F1 + ... + xm*Fm >= 0 
  6. //             Tr F(x) <= M
  7. //
  8. // maximize    -Tr F0*(Z-zI) - Mz
  9. // subject to  Tr Fi*(Z-zI) = c_i
  10. //             Z >= 0, z>= 0
  11. //
  12. // Convergence criteria:
  13. // (1) maxiters is exceeded
  14. // (2) duality gap is less than abstol
  15. // (3) primal and dual objective are both positive and
  16. //     duality gap is less than (reltol * dual objective)
  17. //     or primal and dual objective are both negative and
  18. //     duality gap is less than (reltol * minus the primal objective)
  19. // (4) reltol is negative and
  20. //     primal objective is less than tv or dual objective is greater
  21. //     than tv
  22. //
  23. // Input arguments:
  24. // F:         (sum_i n_i^2) times (m+1) matrix
  25. //            [ F_0^1(:) F_1^1(:) ... F_m^1(:) ]
  26. //            [ F_0^2(:) F_1^2(:) ... F_m^2(:) ]
  27. //                ...      ...          ...
  28. //            [ F_0^L(:) F_1^L(:) ... F_m^L(:) ]
  29. //            F_i^j: jth block of F_i, size n_i times n_i.
  30. // blck_szs:  L-vector [n_1 ... n_L], dimensions of diagonal blocks.
  31. // c:         m-vector.  Specifies primal objective.
  32. // x0:        m-vector.  The primal starting point.  F(x0) > 0.  
  33. // M:         scalar. M > Tr F(x0).    
  34. // nu:        >= 1.0.  Controls the rate of convergence.
  35. // abstol:    absolute tolerance.
  36. // reltol:    relative tolerance.  Has a special meaning when negative.
  37. // tv:        target value.
  38. // maxiters:  maximum number of iterations.
  39. //
  40. // Output arguments:
  41. // x:         m-vector; last primal iterate.
  42. // Z:         last dual iterate; block-diagonal matrix stored as 
  43. //            [ Z^1(:);  Z^2(:); ... ; Z^L(:) ].
  44. // z:         scalar part of last dual iterate.  
  45. // ul:        ul(1): primal objective, ul(1): dual objective.
  46. // iters:     number of iterations taken.
  47.  
  48.  
  49. [rowf,colf]=size(F);
  50. m = colf-1;
  51. if (rowf ~= sum(blck_szs.*blck_szs))
  52.     error('Dimensions of F do not match blck_szs.');
  53. end;
  54.  
  55. [rowx0,colx0]=size(x0);
  56. if (rowx0 ~= m) | (colx0 ~= 1)
  57.    error('x0 must be an m-vector.'); 
  58. end;
  59.  
  60. if (prod(size(x0)) ~= m), 
  61.    error('c must be an m-vector.'); 
  62. end;
  63.  
  64. // I is the identity
  65. I = zeros(rowf,1);
  66. blck_szs=matrix(blck_szs,1,prod(size(blck_szs)));
  67. k=0;  for n=blck_szs,
  68.    I(k+[1:n*n]) = matrix(eye(n,n),n*n,1);   // identity
  69.    k = k+n*n;   // k = sum n_i*n_i
  70. end;
  71.  
  72. // Z0 = projection of I on dual feasible space 
  73.  
  74. Z0 = I-F(:,2:m+1) * ...
  75.      ( (F(:,2:m+1)'*F(:,2:m+1)) \ ( F(:,2:m+1)'*I - c ) );
  76.  
  77. // mineigZ is the smallest eigenvalue of Z0
  78. mineigZ = 0.0;
  79. k=0; for n=blck_szs,
  80.   mineigZ = min(mineigZ, min(real(spec(matrix(Z0(k+[1:n*n]),n,n)))));
  81.   k=k+n*n;
  82. end;
  83.  
  84. // z = max( 1e-5, -1.1*mineigZ )
  85. Z0(k+1) = max( 1e-5, -1.1*mineigZ);  
  86. Z0(1:k) = Z0(1:k) + Z0(k+1)*I; 
  87.  
  88. if (M < I'*F*[1;x0] + 1e-5), 
  89.    error('M must be strictly greater than trace of F(x0).'); 
  90. end;
  91.  
  92. // add scalar block Tr F(x) <= M
  93. F = [F; M-I'*F(:,1),-I'*F(:,2:m+1)]; 
  94. blck_szs = [blck_szs,1];
  95. [x,Z,ul,info]=...
  96.   semidef(x0,pack(Z0),pack(F),blck_szs,c,[nu,abstol,reltol,tv,maxiters]);
  97. iters = info(2);
  98. nz=prod(size(Z))
  99. z=Z(nz)
  100. Z=unpack(Z(1:nz-1),blck_szs(1:prod(size(blck_szs))-1))
  101. Z = Z(1:k);
  102.  
  103.  
  104.  
  105.